#include <Eigen/Dense>
#include <iostream>

Eigen::VectorXd polyfit2(const Eigen::VectorXd& x, const Eigen::VectorXd& y, int degree) {
    // 1. 数据归一化（提升数值稳定性）
    const double x_mean = x.mean();
    const double x_scale = 2.0 / (x.maxCoeff() - x.minCoeff());
    Eigen::VectorXd x_norm = (x.array() - x_mean) * x_scale;

    // 2. 构建范德蒙矩阵（递推计算幂次）
    Eigen::MatrixXd A(x.size(), degree + 1);
    for (int i = 0; i < x.size(); i++) {
        A(i, 0) = 1.0;  // x^0
        for (int j = 1; j <= degree; j++) {
            A(i, j) = A(i, j - 1) * x_norm(i);  // 递推计算x^j
        }
    }

    // 3. 列主元QR分解求解
    Eigen::VectorXd coeffs_norm = A.colPivHouseholderQr().solve(y);

    // 4. 将系数转换回原始尺度
    Eigen::VectorXd coeffs(degree + 1);
    coeffs.setZero();
    for (int j = 0; j <= degree; j++) {
        for (int k = j; k <= degree; k++) {
            double comb = 1.0;  // 组合数C(k,j)
            for (int i = 1; i <= j; i++) comb *= (k - j + i) / double(i);
            coeffs(j) += coeffs_norm(k) * std::pow(-x_mean * x_scale, k - j) *
                std::pow(x_scale, j) * comb;
        }
    }
    return coeffs;
}

int test_main3() 
{
    // 您的无噪声数据
    Eigen::VectorXd x(5), y(5);
    x << 1, 2, 3, 4, 5;
    //y << 3, 5, 8, 9, 12;
    y << 1.1, 3.9, 8.2, 15.3, 24.5;

    Eigen::VectorXd coeffs = polyfit2(x, y, 2);
    std::cout << "拟合系数（升幂排列）: " << coeffs.transpose() << std::endl;
    // 输出: [0.1, 0, 1] 对应 0.1 + 0*x + 1*x²

    // 转换为降幂排列（如需与NumPy对齐）
    std::cout << "降幂排列: [" << coeffs(2) << ", " << coeffs(1) << ", " << coeffs(0) << "]\n";
    return 0;
}